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ABSTRACT 

An  Inverse  design  procedure  for  airfoils,  based  on  hodograph  tech¬ 
niques,  has  been  developed.  For  subcritical  flows  the  subsonic  portion 
of  the  pressure  distribution  is  prescribed;  for  supercritical  flows  the 
stream  function  on  the  sonic  line  is  given  instead  of  the  supercritical 
portion  of  the  pressure  distribution.  In  the  special  variables  we  use, 
the  equation  for  the  stream  function  may  be  solved  iteratively  using  a 
fast  Poisson  solver,  for  the  subsonic  portion  of  the  flow;  for  super¬ 
critical  flows,  the  supersonic  portion  is  calculated  using  a  character¬ 
istic  calculation.  The  results  are  then  mapped  back  to  the  physical 
plane  to  determine  the  airfoil  shape.  Both  subcritical  and  supercritical 
results  are  presented.  They  show  good  agreement  with  the  direct  compu¬ 
tation  of  the  flow  past  the  designed  airfoil. 
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INTRODUCTION 

In  the  last  decade  there  has  been  much  interest  in  transonic  aero¬ 
dynamics.  In  the  case  of  conmercial  aircraft  this  interest  stems  mainly 
from  the  requirement  for  increased  efficiency.  This  requirement  has 
driven  the  operating  conditions  of  aircraft  compressors,  turbines,  and 
propellers,  as  well  as  the  aircraft  themselves,  into  the  transonic  regime. 
Unfortunately  once  local  regions  of  supersonic  flow  occur  shock  waves  are 
likely,  with  the  attendant  wave  drag  and  boundary-layer  separation  losses. 

A  shock-free  flow  would,  of  course,  be  desirable  as  it  avoids  these  losses. 
But  in  the  mid  1950s  Morawetz^  showed  that  shock-free,  two-dimensional, 
irrotational  near-sonic  flows  are  mathematically  isolated  (i.e.,  any  small 
arbitrary  changes  in  the  upstream  flow  or  airfoil  shape  will  lead  to  the 
formation  of  a  shock  wave).  As  a  consequence,  it  was  felt  for  some  time 
that  shock-free  flows  could  not  be  achieved  in  practice.  Subsequent  wind 
tunnel  research  carried  out  by  Pearcey,2  Whitcomb  and  Clark3  and  Spee4 
did,  however,  lead  to  the  development  of  practical  "shock-free"  airfoil 
sections.  Parallel  with  these  efforts,  shock-free  supercritical  airfoil 

flows  have  been  generated  in  the  hodograph  plane  by  Nieuwland,5  Bauer, 

6  7  8 

Garabedian  and  Korn,  Boerstoel  and  Huizlng,  and  Sobieczky.  More 

recently,  the  development  of  sophisticated  numerical  codes  for  the 

analysis  of  transonic  flow  fields  has  led  to  the  design  of  both  airfoils 

and  wings  by  numerical  optimization  (see,  e.g.,  Hicks  and  Henne9). 

The  design  procedure  for  shock-free  supercritical  airfoils  developed 
by  Garabedian  uses  the  method  of  complex  characteristics  and  is  mathe¬ 
matically  elegant.  However,  even  In  Its  user-oriented  version6  practical 
use  of  the  method  is  limited  to  especially  trained  designers.  Neverthe¬ 
less,  some  very  useful  test  cases  for  supercritical  airfoils  and  cascades 
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were  given  which  have  served  well  as  examples  for  various  analysis  codes. 

O 

Another  indirect  approach  was  developed  by  Sobieczky.  Here  the  rheo- 
electric  analogy  to  the  compressible  flow  equations  was  used  to  define  a 
number  of  transonic  airfoil  flows.  The  use  of  an  electric  setup,  of 
course,  did  not  provide  for  economical  use  of  the  method,  but  did  lead  to 
a  greater  understanding  of  this  indirect  approach,  which  had  two  fruitful 
consequences:  1)  By  replacing  the  electric  analog  by  a  digital  computer 
using  a  fast  elliptic  solver,  we  could  use  all  the  experience  of  the  analog 
method  concerning  boundary  conditions  that  result  in  interesting  airfoil 
designs.  2)  With  the  development  of  numerical  algorithms  for  direct  non¬ 
linear  analysis  a  new  method  for  direct  transonic  design  was  developed.^® 
Its  underlying  "fictitious  gas  concept"  is  the  physical  equivalent  of  the 
mathematical  formulation  of  the  indirect  boundary-value  problem.  The 
direct  approach  has  the  advantage  that  it  can  be  used  also  for  three- 
dimensional  wings.  A  description  of  these  two  families  of  methods  applied 
to  airfoils  is  given  In  Ref.  (11). 

The  design  procedure  invoked  here  Is  an  extension  and  improvement  of 
method  1  above;  it  requires  less  than  a  minute  of  CYBER  175  CPU  time  for 
a  supercritical  airfoil,  and  even  less  for  a  subcrltical  airfoil.  While 
the  present  method  is  limited  to  two-dimensional  flow,  it  does  allow  the 
user  to  specify  a  desired  pressure  distribution  and  achieve  It  with  little 
difficulty  through  an  Iterative  process. 

FORMULATION 

The  mathematical  formulation  of  the  problem  assumes  a  steady,  two¬ 


dimensional,  Irrotatlonal  flow  of  a  perfect  gas.  The  basic  equations  of 
motion  are 


t 
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7  •  p£  *  0, 


7  x  £  =  0, 


-4-  =  constant. 


We  satisfy  the  first  of  the  equations  by  introducing  the  usual  stream 


function 


<i»y  =  pu  =  pq  Cos(e), 


♦x  =  -pv  -  -pq  Sin(e) . 


We  satisfy  the  second  by  defining  a  potential  function,  such  that 


v<t>  =  g. 


u  ■  <j>x.  v  *  <py. 


Equations  (4)  and  (5)  can  be  combined  to  give 


<i>  *  —  \i>  , 

*  P  7 


A  * - 6  , 

vy  p  yx 

which  can  be  reduced  to  a  single  equation  In  either  ^  or  <f>,  viz., 

p  P 

*xx  +  ^yy  "  {f}  *x  +  3 1  V 

^xx  +  ^yy  *  "^"p”3  ^x  “  ^ ~p  3  ^y* 


(7) 
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wi  th 

p  «  p(q)»  q  =  |v<j>| .  (8) 

Equations  (7)  are  elliptic  for  M  <  1  (subsonic  flow),  hyperbolic  for  M  >  1 
(supersonic  flow),  and  parabolic  where  M  =  1  (sonic  flow).  Figure  (1)  is 
a  sketch  illustrating  a  shock-free  supercritical  flow  field  depicting 
subsonic  and  supersonic  regions  of  the  flow. 

To  solve  either  of  equations  (7)  two  boundary  conditions  must  be  pro¬ 
vided.  The  first  of  these  represents  the  behavior  of  the  flow  in  the  far 
field;  the  second  represents  the  flow  tangency  condition  on  the  airfoil 
surface.  In  the  inverse-design  procedure  under  study  here,  the  pressure 
distribution  is  used  as  an  input  rather  than  the  airfoil  coordinates  and 
our  goal  is  to  find  the  shape  that  an  airfoil  must  have  to  achieve  this 
input  pressure  distribution.  Accordingly  for  the  inverse  problem,  neither 
a  Neumann  boundary  condition  for  <J>  nor  a  Dirichlet  boundary  condition  for 
if/  can  be  prescribed  on  the  airfoil  surface  since  it  is  unknown.  This  is 
in  contrast  to  direct  analysis  methods  where  the  airfoil  geometry  is  known 
prior  to  computational  procedure.  The  inverse  problem  of  a  prescribed 
pressure  on  a  given  airfoil  Is  ill -posed.  Here  we  prescribe  a  general 
pressure  distribution  and  find  an  airfoil  that  has  a  pressure  distribution 
like  that  prescribed. 


The  Hodograph  Transformation 

In  two-dimensional  Irrotational  flow  the  nonlinear  equations  (7)  for 
steady  flow  can  be  rendered  linear  by  changing  the  role  of  the  dependent 
and  independent  variables.  We  Introduce  the  complex  velocity 


iv  *  qe" 


I 

* 
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with  u,  v,  q  and  6  being  functions  of  a  complex  variable 


z  *  x  +  1y. 


We  then  have,  using  equations  (6,  9)  in  differential  notation. 


+  i  d^  *  (u  -  1v)(dx  +  idy) 
=  qe'i8  dz 


and  hence,  with  q  and  0  as  independent  variables,  and  p  =  p ( | V<p | )  we  find 


(iUL)  =  ® _ (1£  +  4  1  M\ 

vae;q  q  '30  1  p  30;’ 


(M)  =  e _ (M  +  i  1  li) 

v3q;0  q  '3q  p  3q;* 


Differentiating  the  first  equation  of  (12)  with  respect  to  q  and  the  second 
with  respect  to  0  and  then  equating  real  and  imaginary  parts,  we  obtain  the 
hodograph  equations,  viz.. 


li ,  i 

30  p  3q  ’ 


li.  =  _  _L  (1 
3q  pq  U 


-  M2!  II 


It  is  essential  in  the  present  design  procedure  to  introduce  the 
Prandtl -Meyer  function,  v,  defined  by 


L/irn?Tf 


in  place  of  the  velocity  q  in  equations  (13),  which  then  take  their  canonical 
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form  and  become 


*Q  =  K(v)  \|>v, 

<PV  s  ±K(v)  . 


(15) 


Here  the  ±  signs  refer  to  supersonic  and  subsonic  conditions  respectively  and 
K(v)  •  K[v(M)]  =  [|1  -  M2|]1/2/p[q(M)].  (15a) 


A  typical  hodograph  representation  (with  v,  0  as  independent  variables)  of 
the  flow  field  sketched  in  Figure  (1)  is  depicted  in  Figure  (2a). 

The  airfoil  maps  into  a  closed  curve  containing  the  stagnation  point  S 
at  infinity  and  the  sonic  line  at  points  a,  b.  The  analytical  structure  of 
ip  near  v  =  0  and  a  smooth  curvature  of  the  airfoil  at  the  sonic  line  require 
the  local  structure  of  the  airfoil  contour  near  the  0-axis  as  sketched  in 
Figure  (2b).  This  weakly  singular  behavior  is  of  importance  in  the  solution 
of  the  boundary-value  problem,  with  ip  given  along  the  elliptic  boundary 
(shaded  line  in  detail  of  Figure  2).  The  region  of  the  flow  in  Figure  (1) 
bounded  by  the  dotted  contour  and  the  airfoil  represents  a  region  in  which 
every  point  has  a  velocity  and  flow  angle  equal  to  that  of  some  other  point 
outside  this  region,  e.g.,  qg  =  q^  and  0g  =  e^.  Therefore  points  e  and  f 
will  correspond  to  the  same  point  in  the  hodograph  plane,  which  must  be 
considered  as  a  Riemann  surface  of  two  sheets  with  a  branch  cut  (lines  dn, 

1 7 

cn)  connecting  them.  The  local  nature  of  this  flow  was  studied  by  Lighthill 
In  his  hodograph  study  of  compressible  flows  around  lifting  airfoils. 

Because  the  governing  equations  (15)  are  linear  in  the  hodograph  plane, 
there  Is  usually  no  particular  difficulty  in  finding  solutions  to  them,  by 
numerical  methods  if  necessary.  However,  the  presence  of  the  second  sheet 
of  the  two-sheeted  Riemann  surface  and  the  apriori  unknown  location  of  the 
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airfoil  surface  represent  major  obstacles  in  solving  the  governing  equations 
in  this  plane.  From  equations  (15)  it  can  be  shown  that  the  second-order 
derivatives  for  both  and  <j>  either  are  the  Laplacian  or  the  wave  operator 
depending  on  whether  or  not  the  flow  is  subsonic  or  supersonic.  Thus  the 
equations  for  the  subsonic  flow  are  invariant  in  form  under  a  conformal 
transformation. 


Subsonic  Flow  Domain 

We  proceed  by  assuming  a  conformal  map  of  the  subsonic  portion  of  the 
two-sheeted  Riemann  surface  of  Figure  (2)  into  the  unit  circle  of  Figure  (3). 
(The  double-connected  infinite  domain  is  mapped  into  a  finite  simply  con¬ 
nected  domain  through  an  exponential  mapping  followed  by  a  square  root 
mapping  to  unfold  the  Riemann  sheets.)  Here  part  of  the  boundary  of  the 
unit  circle  corresponds  to  the  airfoil  surface  wetted  by  subsonic  flow;  the 
remaining  part  corresponds  to  the  sonic  line.  The  segment  comprising  the 
sonic  line,  <  u>  <  u>a  is  chosen  and  the  Mach  number  (or  pressure)  on  the 
subsonic  part  of  the  airfoil  prescribed.  On  the  sonic  line  segment  the 
Mach  number  is  of  course  1.  8ecause  Bernoulli's  equation,  viz.. 


+ 


constant 


provides  a  correspondence  between  the  values  of  p  and  the  local  sound 
speed,  we  may  formulate  the  inverse-design  problem  in  terms  of  either  M  or 
p,  and  our  choice  of  M  is  only  a  matter  of  convenience.  A  typical  choice 
of  M  for  a  subcritical  airfoil  with  a  cusped  trailing  edge  is  illustrated 
in  Figure  (4). 

With  the  Mach  number  given  on  the  boundary  of  the  unit  circle,  and 
with  the  subsonic  portion  of  the  flow  inside  the  circle,  we  take  advantage 
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of  the  fact  that  the  mapping  to  the  CQ-plane  is  conformal.  Thus,  the 
Prandtl-Meyer  function  v,  and  the  flow  deflection  angle,  0,  are  conjugate 
harmonics,  i.e., 

F(£0)  =  v  +  10 
or 


V2v(5q)  =  0, 

(16) 

v20(co)  =  o. 

(17) 

Here  F(£Q)  is  the  mapping  function  and  =  re1w  with  r  and  w  being  the 
radial  and  angular  coordinate  measured  in  the  £Q-plane.  Boundary  conditions 
for  equation  (16)  are  provided  through  the  use  of  equation  (14)  relating  v 
and  M;  that  is,  knowing  the  Mach  number  distribution  on  the  unit  circle,  we 
calculate  v  employing  equation  (14).  We  then  solve  Laplace's  equation  for 
v  inside  the  unit  circle  using  Fourier  series,  which  accordingly  determines 
the  flow  deflection  angle  0  to  within  a  constant.  However,  the  Prandtl- 
Meyer  function  v  is  logarithmically  singular  in  q(v  «  log(q))  at  the 
stagnation  point,  S,  which,  for  convenience,  is  positioned  at  £  •  -1  in 
the  SQ-plane.  Therefore,  In  order  to  solve  the  boundary  value  problem  for 
v  with  Fourier  series  we  need  first  to  subtract  the  logarithmic  behavior 
at  point  S.  This  Is  done  as  follows: 

let  F(?0)  -  FU0)  -  log«-0  +  1), 

G (50)  *  v(C0)  -  Re{log(C0  +  1)>, 

H(£0)  ■  6(S0)  -  Im{logU0  +  D), 


then  equations  (16,  17)  become 


10 


v2g(c0)  =  o. 

(18) 

V2H(?0)  =  0. 

(19) 

Here  Re{ - }  and  Im{ - }  are  the  real  and  imaginary  parts  of  log(£Q  +  1) 

respectively.  Equation  (18)  is  then  solved  inside  the  unit  circle  using 
Fourier  series  subject  to  the  following  boundary  condition 

G(r=l ,w)  =  v(l,u>)  -  Re{log(eia)  +  1)}. 

Having  obtained  the  solution  for  G(r,(ij)  in  the  unit  circle,  we  then  add  back 
the  logarithmic  singularity  to  preserve  the  singular  behavior  of  v  at  the 
stagnation  point. 

Reformulating  the  definitions  for  the  partial  derivatives  in  equations 
(15)  in  terms  of  r  and  w  we  obtain 

‘V  *  '  r  K(v)  V 

and  (20) 

<(>u  =  r  K(v)  i pr. 

Equations  (20)  can  be  reduced  to  a  single  equation  in  either  ^  or  $. 
Eliminating  <j>  through  cross  differentiation  of  equations  (20),  we  find  the 
equation  for  if* 

r\r  *  %  +  *  f(M){r2v^r  +  vjj.  (21) 

Here  K  is  a  function  of  M  through  equation  (15a),  M  is  a  function  of  v 
through  equation  (14),  and 


.2 1-3/2 


(21a) 


n 


Equation  (21)  is  the  CQ-p1ane  counterpart  of  the  physical -plane  equation 
(7)  for  the  stream  function.  The  transformation  to  hodograph  variables  v, 

6,  followed  by  a  conformal  transformation  to  the  £Q-plane  results  in  a  linear 
second-order  partial  differential  equation  for  r|i.  In  addition  to  this 
linearity,  the  major  advantages  of  these  two  transformations  are  the  unfold¬ 
ing  of  the  two-sheeted  hodograph  surface  to  a  single  sheet,  and  the  repre¬ 
sentation  of  the  subsonic  portion  of  an  unknown  airfoil  by  a  unit  circle.  But  an 
these  advances  are  not  without  attendant  complexities,  albeit  minor  ones. 

These  difficulties  include  the  presence  of  point  "I"  Figure  (3)  representing 
the  far  field  inside  the  unit  circle,  at  which  the  stream  function  is 
singular,  and  the  singularity  in  f(M)  in  equation  (21a)  at  H  *  1.  The  first 
of  these  minor  difficulties  is  circumvented  through  a  coordinate  transforma¬ 
tion  of  the  £;Q-plane  into  the  5-plane  sketched  in  Figure  (5).  The 
transformation  employed  is  the  bilinear  transformation  which  has  the  form 


Here 


5  =  (i-_ S-)  e-1'8 

*  6  +  Of  '  e 


A  «  (L  -  1), 


■°I  °I 


and 


B  ’  5o,(1  -  50l>’ 


C=  !50  -  1). 


o  ■  o  -  e0i>. 


(22) 


The  second  difficulty  is  the  singular  behavior  at  f(M)  at  M  *  1,  which 
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-1/3 

implies  that  if>r  *  v  '  on  the  sonic  line.  This  is  circumvented  by  assuming 
that  the  Mach  number  on  the  sonic  line  is  slightly  less  than  1,  e.g., 

Msonic  =  0*995  and  using  the  local  behavior  of  \p  to  extend  the  results  to 
M  »  1 . 

With  r  and  w  as  new  independent  variables,  the  governing  equation  for 
1 1»  becomes 


f2*ff 


+  ♦(jj;  ■  f(M){r 


0)  w 


(23) 


In  equation  (22)  the  quantities  v*  and  are  evaluated  through  the  utiliza¬ 
tion  of  the  solution  obtained  for  equation  (16)  in  conjunction  with  equation 
(22).  That  is  to  say 


V'  -  7  Re(CFU)). 

r 

Vjj  *  -Im(CFU)}, 

where 


F(5)  -  e”1u(v  -  i  v  ) 

r  r  W  (B  -  DE  e,p) 


BC  +  DA  _iB 
- TITe  • 


On  the  portion  of  the  airfoil  surface  wetted  by  subsonic  flow  we  have 


<Hl  iu)  *  0  d>a  >.  to 

On  the  remaining  portion  of  the  unit  circle,  i.e.,  the  sonic  line,  an 
arbitrary  distribution  for  the  stream  function,  is  prescribed,  viz.. 


<p(l  ,ui)  *  |J>S  a  g(u>)  U)a  <  u>  <  d)^ 


The  compressible  far-fleld  stream  function,  ,  is  related  to  its  incom- 

00 

presslble  counterpart,  i'  ,  via 

00 


(24) 


1  q»e1Ct  ir 

=  <T(ry  Im{(c0  -  s  )  "  2?  log(co  '  s0l)}- 

Equation  (24)  is  a  consequence  of  the  fact  that,  at  point  "I"  which  presents 
infinity  in  the  5Q-plane  equations  (20)  reduce  to  the  familiar  Cauchy- 
Riemann  equations  (or  incompressible  flow  equations)  expressed  in  polar 
coordinates.  Thus  the  complex  potential  defining  the  far-field  stream 
function  is  that  for  an  incompressible  flow,  ¥  ,  which  is  related  to  its 
compressible  counterpart,  ij^,  via  equation  (24).  Equation  (24),  in  conjunc¬ 
tion  with  equation  (22),  is  then  employed  in  expressing  the  far-field 
boundary  condition  in  the  £-plane. 

The  boundary  value  problem  for  the  stream  function  is  now  complete,  and 
equation  (22)  is  solved  iteratively  using  a  sixth-order  accurate  fast 

1 3 

Poisson  solver  devised  by  Roache.  The  stagnation  streamline  leaves  the 
contour  of  an  airfoil  at  a  cusped  or  wedged  trailing  edge.  This  defines,  in 
direct  analysis,  the  amount  of  circulation  around  the  airfoil  and  the  loca¬ 
tion  of  the  stagnation  point  near  the  leading  edge.  For  an  indirect  method 
such  as  the  present  one,  the  situation  is  exactly  the  reverse:  Since  the 
(mapped)  stagnation  point  location  is  given,  we  have  to  vary  the  circulation 
until  the  stagnation  streamline  reaches  the  surface  exactly  at  the  stagna¬ 
tion  point.  The  mapped  location  T  of  the  trailing  edge  results  from  this 
adjustment  and  the  trailing  edge  shape  Is  then  cusped:  In  the  computational 
plane  the  stagnation  streamline  leaves  the  airfoil  contour  at  T  at  a  90° 
angle. 

The  results  for  the  gradients  and  ^  along  the  elliptic  boundary 
consisting  of  the  subsonic  portion  of  the  airfoil  and  the  sonic  line  are 


then  utilized  in  finding  the  inverse-map  to  the  physical  plane,  equation 
(11),  which  can  be  written  as 
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dx(u>)  =  ^  {K( v)Cos(0 )\ps  -  ^  Sin(9)^}dc3, 
and  (25) 

oy(&)  s  ^  {K(v)Sin(0)ip<;  +  1  Cos(0)^}d(j3. 


Equations  (25)  are  numerically  integrated  along  the  unit  circle  starting 
from  the  trailing  edge  T  on  the  airfoil's  lower  surface  moving  toward  the 
stagnation  point  S  and  finally  ending  at  the  trailing  edge  T  on  the  air¬ 
foil's  upper  surface,  i.e.. 


(26) 


The  resulting  airfoil  configuration  is  then  checked  to  see  if  it  has  a 
reasonable  thickness  distribution.  If  not,  then  the  input  design  parameters 
must  be  altered  and  the  design  procedure  repeated.  When  a  suitable  airfoil 
configuration  is  found,  the  second  of  equations  (20)  is  then  numerically 
integrated  to  find  the  potential  distribution  along  the  elliptic  boundary, 
i.e. , 


2ir +uly 

4>(1»<*>)  a  <J>0O.Wj-)  +  |  K(v)i^  da!. 


(27) 
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Results  obtained  from  the  solution  of  the  subsonic  flow  equations  (20) 
on  the  sonic  line  are  then  used  as  initial  values  for  solving  the  supersonic 
flow  equations  (15)  in  the  hodograph  plane  using  the  method  of  character¬ 
istics  explained  below. 


Supersonic  Flow  Domain 

It  follows  from  equation  (15)  for  the  flow  in  the  embedded  supersonic 
region,  where  M  >  1,  that  both  ip  and  $  satisfy  the  linear  wave  equations 

*vv  -  *ee  =  -  kU  V 

and 

♦w  ‘  *06  *  _K(v)  K(W]  V 

which  have  characteristics  of  the  form 

&**'•  <28> 
Introducing  the  characteristics  coordinates  C,  n  defined  by 


5  *  v  +  0, 


(28a) 


n  =  u  -  6, 


and  reformulating  the  partial  derivatives  In  equation  (15)  in  terms  of  the 
new  coordinates,  we  obtain 


<f>£  *  K(v) 

and  (29) 

♦n  *  -K(v)  *n. 

Equations  (29)  may  also  be  expressed  as 
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on  C  *  constant  d<j>  +  K(v)<fy  =  0, 


and  (30) 

on  n  *  constant  d<j»  -  K(v)d4>  =  0. 


With  the  stream  functions  ^$(0)  and  the  potential  <t>s(e)  data  known  on 
the  sonic  line  with  coordinates  xg ( 9)  and  ys(e),  we  proceed  to  solve  the 
compatability  equations  (30)  holding  along  the  two  families  of  characteris¬ 
tics  using  a  step-by-step  numerical  scheme  (Massau  finite  difference  scheme). 
The  basic  concept  of  this  scheme  is  the  following: 

Through  each  point  A,B,C,D,E,F,  of  the  sonic  line  (line  AF  in  Figure 
(6))  two  Mach  lines  pass,  one  of  the  first  family,  n  *  constant  (a^,  b^ ,  c^ , 

. ..fj),  and  one  of  the  second  family,  £  =  constant  (a2,  b2>  c2,...f2). 

Since  the  flow  properties  at  the  points  of  the  line  AF  are  known,  the  con¬ 
stants  of  equations  (30)  for  each  characteristic  line  a-j...^,  a2...f2,  at 
each  point  of  AF,  are  also  known.  Consequently,  equations  (30)  applied 
along  two  characteristics  of  opposite  families,  e.g.,  b-j  and  a2  of  Figure 
(6),  give  two  equations  relating  <j>Q  to  These  can  be  solved  to  obtain 
<j>  and  \p  at  point  G  as  functions  of  $  and  ^  at  the  points  A  and  B,  viz.. 


^G  *  2  **A  +  V  +  2KTvT  (*A  ’  V* 


(31) 


*G  "  \  +  V  +  **A  "  V* 


In  (31),  If  we  think  of  v  as  a  specified  function  of  £  and  n  on  a  rectangu¬ 
lar  grid,  the  coefficient  K(v)  may  be  approximated  by  average  values  such 
that  the  numerical  step-by-step  scheme  becomes  second-order  accurate  in  the 
mesh  size,  h.  The  flow  deflection  angle,  0,  and  the  Prandtl -Meyer  function, 
v,  are  evaluated  at  the  grid  nodes  using  (28a). 
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Having  solved  for  the  flow  properties  at  all  grid  nodes,  we  then 
proceed  to  search  for  points  of  zero  stream  line  values  as  they  represent 
points  on  the  airfoil  upper  surface. 

The  inverse  transformation  from  the  hodograph  plane  to  the  physical 
plane  is  given  by  equation  (11),  which  we  can  write  at  every  grid  point  as 

dx  =  £21121  4 . a*, 

q  pq 

(32) 

dy  .  4  t  COJM  d„. 

The  Jacobian  of  transformation  is 


“  ~T  (V0  '  Vv5'  (33) 

pq 

Utilizing  equations  (15),  we  find  that  for  supersonic  flow 


which  indicates  the  possibility  of  a  vanishing  Jacobian  whenever  |ij>J  * 

It  can  be  easily  shown  that  this  condition  is  satisfied  whenever  a 
•  constant  line  is  tangent  to  the  characteristics  in  the  hodograph  v-e 
plane,  i.e.,  whenever 


v° 

or 

*  o. 

The  locus  of  points  in  the  hodograph  plane  for  which  J  -  0  is  known  as  a 
limit  line.  In  the  physical  plane  the  image  of  such  a  line  is  a  cusped 
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curve  along  which  the  physical  surface  can  be  thought  to  be  folded  upon  it¬ 
self.  Such  limit  lines  indicate  the  inconsistency  of  the  sonic  line  data 
with  shock-free  flow.  Thus,  they  are  only  acceptable  if  they  occur  below 
the  airfoil's  upper  surface.  However,  if  they  occur  above  or  on  the  air¬ 
foil's  upper  surface,  the  sonic  line  data  must  be  altered  by  adjusting  the 
input  subcritical  pressure  or  the  sonic  line  stream  function  distributions 
or  both,  and  the  design  procedure  repeated. 

The  Use  of  the  Computer  Program  Package 

The  computer  program  package  consists  of  sixteen  FORTRAN  IV  sub¬ 
routines  performing  calculations  and  data  transfer.  The  package  makes  use 
of  two  permanent  files  for  results  to  be  retained  and  one  file  for  input 
data. 

For  a  specific  airfoil  design  the  program  has  to  be  operated  in  the 
following  manner: 

a  —  Using  the  input  Mach  number  and,  for  supercritical  airfoils,  the  sonic 
line  stream  function  distributions,  both  of  which  must  be  reasonable, 
in  addition  to  the  two  design  parameters  M^  and  a,  we  obtain  a 
subsonic  boundary  configuration  and  have  it  displayed  on  a  computer 
terminal  with  graphic  capabilities. 

b  —  The  level  of  the  input  Mach  number  distribution  on  the  lower  surface 
(0  £  a)  <  ir)  is  moved  slightly  In  order  to  shift  the  location  of  the 
far-field  singularity  In  the  computational  plane.  This  in  turn  is 
sufficient  to  reduce  the  x-gap  of  the  trailing  edge  of  the  resulting 
configuration  to  a  preset  value  (e.g.,  10  ). 
c  —  The  Input  Mach  number  distribution  Is  then  altered  in  the  vicinity  of 
the  stagnation  point  of  the  resulting  configuration  to  reduce  the 
y-gap  at  that  point. 
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d  —  Once  an  acceptable  subsonic  boundary  configuration  (displayed  on  the 
terminal)  is  obtained  the  supersonic  flow  field  computation  is  then 
carried  on  to  determine  the  airfoil  coordinates  under  the  sonic  line. 

If  a  limit  line  is  encountered  in  the  supersonic  flow  domain,  steps 
a  through  d  are  repeated  for  modified  input  Mach  number  and  sonic 
line  stream  function  distributions. 

Figure  (7)  is  a  flow  chart  illustrating  the  design  procedure  discussed  above. 

COMPUTATIONAL  RESULTS 

In  this  section,  a  number  of  numerical  computational  results  are  pre¬ 
sented,  and  some  experience  with  the  algorithm  is  discussed.  A  few  examples 
of  airfoil  designs  are  presented  in  Figures  (8)  through  (20).  The  airfoils 
cover  a  range  of  free-stream  Mach  numbers  and  lift  coefficients.  A  comparison 
with  results  obtained  from  the  direct  computation  of  the  flow  field  utilizing 
the  designed  airfoil  geometry  as  input  is  also  illustrated  in  Figures  (16) 

and  (20)  for  a  subcritical  and  a  supercritical  airfoil  respectively.  The 

14 

analysis  code,  FL06,  devised  by  Jameson  utilizes  a  finite-difference 
method  to  solve  the  governing  equations  iteratively  after  mapping  the  airfoil 
into  the  unit  circle. 

For  the  design  of  subcritical  airfoils  with  cusped  trailing  edges. 

Figures  (8)  through  (15)  illustrate  the  different  input  Mach  number  distribu¬ 
tions  required  for  the  design  procedure  and  the  resulting  airfoil  shapes  at 
M^  *  0.573,  0.589,  and  0.622.  Through  a  study  of  the  above  airfoils  we  have 
demonstrated  that  small  modifications  in  the  input  Mach  number  distribution 
near  the  stagnation  point  have  a  substantial  influence  on  the  y-gap  at  the 
airfoil  trailing  edge  and,  hence,  on  airfoil  closure.  That  is.  Increasing 
the  slope  of  the  Mach  number  distribution  at  the  stagnation  point  results 
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in  an  airfoil  with  a  smaller  nose  radius  and  consequently  a  smaller  thickness 
to  chord  ratio;  this  ultimately  results  in  a  smaller  y-gap  at  the  airfoil 
trailing  edge.  Experience  has  shown  that  in  the  process  of  minimizing  the 
y-gap  at  the  trailing  edge  the  level  of  the  input  Mach  number  distribution 
has  to  be  altered  on  the  lower  surface  or  the  upper  surface,  or  both,  to 

_3 

maintain  the  x-gap  at  a  present  value,  e.g.,  10  . 

Figure  (16)  compares  the  results  obtained  from  analysis  computation 
of  the  flow  field  utilizing  the  designed  airfoil  geometry  shown  in  Figure  (9) 
as  input.  Results  obtained  from  both  the  design  and  analysis  show  good 
agreement  except  near  the  trailing  edge. 

The  input  Mach  number  distribution  for  a  slightly  supercritical  airfoil 
with  a  cusped  trailing  edge  at  M^  =  0.642  is  shown  in  Figure  (17).  Results 
obtained  from  the  two  stages  of  the  design  procedure,  i.e.,  computation  of 
the  subsonic  boundary  configuration  and  the  complete  airfoil  shape,  are 
illustrated  in  Figures (18,19) .  It  should  be  mentioned  here  that  while  the 
y-gap  at  the  trailing  edge  for  subcritical  airfoils  can  be  controlled  easily, 
more  effort  is  required  in  the  case  of  supercritical  airfoils.  This  is  due 
to  the  representation  of  the  sonic  line  as  part  of  the  subsonic  boundary 
configuration  (see  Figure  (18)). 

Figure  (20)  illustrates  a  comparison  of  results  obtained  from  the  direct 
computation  of  the  flow  field  of  a  supercritical  airfoil  with  cusped  trailing 
edge  at  M^  »  0.75  resulting  from  the  design  procedure  Invoked  here  with  a 
different  solver  for  the  stream  function  equation.  It  can  be  noticed  that 
good  agreement  is  obtained  for  the  subcritical  pressures.  The  slight 
disagreement  between  the  supercritical  pressures  that  Is  observed  is  due 
to  the  numerical  difficulty  in  computing  the  normal  derivative  of  the  stream 


function  along  the  sonic  line;  this  has  a  direct  effect  on  the  sonic  line 
coordinates  and  the  supercritical  pressures. 

CONCLUSION 

Based  on  hodograph  theory  an  efficient  design  procedure  has  been 
developed  that  can  be  used  for  the  aerodynamic  design  of  subcritical  as 
well  as  shock-free  supercritical  airfoil  sections.  An  analysis  of  the 
numerical  solutions  obtained  with  this  design  procedure  indicates  that  only 
one  airfoil  configuration  is  obtained  for  a  desired  input  target  pressure 
and  free-stream  flow  conditions.  Experience  has  also  shown  that  the  input 
Mach  number,  the  free-stream  Mach  number  and  the  sonic  line  stream  function 
distribution  cannot  be  prescribed  independently  of  each  other  if  a  practical 
airfoil  configuration  is  sought.  Subcritical  airfoils  with  a  target  pressure 
distribution  and  shock-free  supercritical  airfoils  with  a  target  subsonic 
pressure  distribution  can  be  found  with  modest  effort  when  the  algorithm 
is  in  the  hands  of  an  experienced  user. 

Within  the  framework  of  our  design  procedure  it  seems  possible  that  the 
procedure  can  be  extended  further  to  study  a  large  class  of  aerodynamically 
interesting  airfoils  and  their  corresponding  pressure  distributions.  Moreover, 
while  the  results  are  given  for  inviscid  flow,  the  same  procedure  can  be 
employed,  iteratively,  with  a  boundary-layer  calculation  (utilizing  airfoils 
with  an  open  trailing  edge)  in  order  to  achieve  viscous  airfoil  designs.  For 
supercritical  airfoils,  the  shock-free  pressure  field  should  make  the  boundary- 
layer  calculation  reliable. 

From  the  engineering  point  of  view,  the  method  was  found  to  be  sufficiently 
accurate  and  flexible  providing  results  at  moderate  cost.  Indeed,  less  than 


a  minute  of  CYBER  175  CPU  time  is  required  for  supercritical  airfoils; 
even  less  is  required  for  subcritical  airfoils. 
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Figure  1.  Sketch  of  shock-free  flow  past  a  lifting 
airfoil  (stagnation  point  S,  trailing 
edge  T,  far-field  I,  branch  point  n). 


Figure  2a 


Figure  2b 


Figure  2.  Hodograph  representation  of  the  flow 
field  sketched  in  figure  (l). 
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Figure  4,  Typical  Input  Mach  number  distribution. 

The  location  of  the  trailing  edge  T  is 
not  known  in  advance. 


Figure  6.  Computation  of  the  supersonic  flow  field  in 
the  hodogreph  plane  (v-9)  using  the  method 
of  characteristics. 
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Figure  9.  Airfoil  818657 
M.  *0,5730  ^  -0.6932  Cjj  -0.0003  A.O 


INVERSE  DESIGN  ALGORITHM 
INPUT  MACH  NUMBER  DISTRIBUTION 
Airfoil  818457 


Angular  position  in  circle  plane  (*•>/  2 *  radians) 


Figure  10,  Input  Mach  number  distribution  for  a  sub- 
critical  airfoil  with  cusped  trailing  edge 
at  Km  -0.573. 


Angular  position  in  circle  plane  (<"/  2* radians) 


Figure  12.  Input  Mach  number  distribution  for  a  sub- 
critical  airfoil  with  cusped  trailing  edge 
at  M—  *0.589. 
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INVERSE  DESIGN  ALGORITHM 
INPUT  MACH  NUMBER  DISTRIBUTION 
Airfoil  809762 
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Figure  14.  Input  Mach  number  distribution  for  a  sub- 
critical  airfoil  vith  cusped  trailing  edge 
at  M^  *0,622. 
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Figure  15.  Airfoil  809762 
M..  -0.622  ■0.7984  CD  -0,0001  A.O.A.  -1.027' 


Mach  Number  (  M  ) 

0.00  0.25  0.50  0.75  1-00  1.25  1.50 


Figure  1?.  Input  Mach  number  distribution  for  a  super¬ 
critical  airfoil  with  cusped  trailing  edge 
at  nm  -0.642, 
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Figure  18.  Airfoil  811064 
M*-  0.642  CL  -  0.943  - 
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Figure  19.  Airfoil  811064 
M„-  0.642  CL  -  0.987  C0  -  -0.0009  A.O.A.  - 
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